A hybrid model for Rydberg gases including exact two-body correlations 
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A model for the simulation of laser-driven interacting Rydberg gases is discussed. Our hybrid 
approach combines an exact two-body treatment of nearby pairs of atoms with an effective treatment 
for separated pairs. A time-independent Monte Carlo technique is used to efficiently determine the 
steady state of the system. The hybrid model predicts features in the pair correlation function arising 
from multi-atom processes which existing models can only partially reproduce. Our interpretation 
of these features shows that higher-order correlations are relevant already at low densities. 
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The distinctive properties of Rydberg atoms [l|, Q ren- 
der them a powerful implementation of a tunable strongly 
interacting quantum many-body system. The signature 
effect of Rydberg atoms is the dipole blockade [3|-[8|, 
which can be understood using rather basic theoretical 
models But ongoing experimental progress pro- 

vides means to explore the correlations on a deeper level. 
Accordingly, more involved observables such as the Man- 
del Q parameter 111, HI |, the pair correlation function [l3| 
or quantum optical effects in the presence of Rydberg in- 
teractions 14, 1^ moved into the focus of interest. State- 
of-the-art experiments push existing models for Rydberg 
gases to the limits of their validity ranges [13, [Tel, and re- 
quire exceedingly long simulation times. Therefore, bet- 
ter simulation techniques are highly desirable. 

One approach to simulate the many-body system is to 
truncate the otherwise exponentially growing state space, 
e.g., based on the Rydberg blockade. Such calculations of 
the exact Hamiltonian dynamics in a truncated Hilbert 
space provide valuable insights [ol, [l7l - [l9| . But they are 
not well suited to model weakly interacting or low-density 
gases, as then the state space becomes too large. Due to 
the treatment on the level of the wave function, incoher- 
ent effects like spontaneous emission or dephasing cannot 
be incorporated in a straightforward way. Also, two-step 
excitation schemes cannot be described, as the interme- 
diate state would spoil the state space reduction. 

An alternative approach is to approximate the inter- 
atomic correlations. The simplest approach is the mean 
field approximation (AiBj) ^ {Ai){Bj) for operators 
A, B acting on two different atoms i, 7, which already 
describes the Rydberg blockade well [sl, [lO^. However, it 
fails for higher correlations and for more involved observ- 
ables [l3, [l5| . A cluster expansion to higher-order atomic 
correlations could successfully describe an experiment on 
coherent population trapping in Rydberg atoms [l^ , but 
this method turned out to give inconsistent results at 
high densities Alternatively, a rate equation model 
was introduced in [2lj. Using an adiabatic elimination 
of the atomic coherences, an effective rate equation for 
the atomic populations alone could be derived. The re- 



sulting equation system still grows exponentially in the 
number of atoms, but can be solved using Monte Carlo 
techniques. This method is capable of modelling a multi- 
step excitation [20]. However, due to the simplified treat- 
ment of the interaction, only single-atom processes can 
be described and correlations are approximated. 

In this work, we discuss a simulation technique combin- 
ing higher predictive power with fast calculation times. 
In our model, atom pairs with distances below a char- 
acteristic length scale Lcr are treated exactly, resulting 
in accurate two-body correlations. Atoms with distances 
above Lcr are incorporated via effective detunings. We 
show that this hybrid approach describes three-body sys- 
tems well over the whole range of interaction strengths. 
To efficiently determine the many-body steady state, we 
eliminate the time-consuming calculations of transition 
rates in each step based on a time-independent Monte 
Carlo technique. We show the consistency of our ap- 
proach with the relevant existing models for simple ob- 
servables as the excitation fraction of a disordered Ry- 
dberg ensemble. In contrast, the hybrid model predicts 
structures in the pair correlation function, which the ex- 
isting models cannot fully reproduce. These structures 
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FIG. 1. (Color online) (a) Our model divides atoms into pairs 
and single atoms allowing an exact treatment of the two-body 
interaction up to a inter-atomic separation Lcr- Encircled 
subsets interact with each other via an effective detuning, 
(b) On a lattice structure, overlapping pairs have to be used. 
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FIG. 2. (Color online) Three models consisting of lower order 
subsets are investigated as simplified models for the = 3 
case. The relative deviation - pir''|/pir'' from the 

exact solution is shown for varying detuning A and interaction 
V with the third atom. Models including the exact treatment 
of pairs give more precise results than the rate equation. The 
dashed line depicts the models which are suited best for a 
given V. Parameters: Q12 - 3 MHz, Q23 - 2 MHz, 721 - 
6 MHz, 732 25 kHz, V12 = 2 MHz. 



arise from higher-order processes, which we show to be 
relevant already at low densities. 

We consider a frozen gas [2] of three-level atoms driven 
by a two-step excitation scheme as shown in Fig. [H 
The lower transition between ground state |1) and in- 
termediate state 1 2) is driven resonantly with Rabi fre- 
quency r^i2, and the upper transition between interme- 
diate state 1 2) and Rydberg state |3) is driven with de- 
tuning A = — uj^2 and Rabi frequency 1^23- 

Before we start with the theoretical analysis, we first 
illustrate the main idea. We start from a two-atom mas- 
ter equation which yields exact results for two atoms. 
But there is no unique method to extend this ansatz 
for pairs to many atoms without using larger atom sub- 
sets. To illustrate this, we consider a system of three 
atoms in which the first two atoms interact strongly with 
coupling V12. A third atom is moved from large dis- 
tance iy := Vi3 = 1^23 = 0) towards the pair until 
the three atoms form an equilateral triangle. For the 
resulting interaction energies V ^ we compare the exact 
three-body steady state solution with different simplified 
models consisting of single atoms or pairs. The result 
for different detunings A is shown in Fig. [2j In the first 
case (a), we consider the single-atom rate equation model 
(SARE) [20I, in which the interactions are absorbed into 
effective detunings for each atom individually, and no 
pairs are formed. The second model (b) accounts for 
the three possible pairs in the calculation. In addition, 
each atom in a pair receives an effective detuning due to 
the interaction with the respective third atom. In the 
third model we assume that V is smaller than F12, and 
treat the first two atoms as pair, and the third as in- 
dividual atom. Again, the mutual interactions between 



single atom and pair are included in effective detunings. 
In Fig. [2] it can be seen that model (a) significantly devi- 
ates from the exact treatment at all interaction strengths. 
Model (b) performs poor at small V ^ because the over- 
lapping pairs mix the exact and approximate interaction 
between two atoms. However, (b) is best at high interac- 
tion strength. The third ansatz (c) reproduces the exact 
steady state in the limit of V ^ by construction, but 
fails at large V, as then correlation with atom 3 cannot 
be neglected. Based on this observation, we construct 
a model in which the two approaches (b) and (c) are 
combined, depending on the mutual interactions of the 
atoms. For this, we introduce a critical interaction (or 
equivalent: a critical distance) that defines when which 
model has to be used. This is indicated by the dashed 
red arrow in Fig. [2l The resulting combination of two 
approaches provides best results over the whole inter- 
action strength range. The increased complexity of our 
model due to the inclusion of exact two-atom correlations 
makes the computation based on established techniques 
impractically slow. We overcome this by a new method 
discussed below. Since this approach involves a combi- 
nation of single-atom and pair descriptions, we call it the 
hybrid model (HM). 

We now proceed with the formal analysis. For a given 
ensemble of N atoms, we introduce a critical distance 
Lcr and use it to divide the set of atoms into pairs 
and single atoms, as schematically shown in Fig. [H Lcr 
should be chosen such that many pairs are included to 
better incorporate two-atom correlations. But on the 
other hand, overlapping pairs should be avoided as seen 
in the = 3 case in Fig. [2l Therefore, we use a 
value which is on the order of the most probable nearest- 
neighbor distance. The single-particle Hamiltonian in 
rotating wave approximation reads 

= l^i2(5« + S^l) + 1^23(4^ + - ^Sfl . (1) 

Here S^^l denotes the operator \a)ii{h\ for the iih atom 
and the Rabi frequencies are assumed to be real. Our 
hybrid model is based on the full two-body Hamiltonian 
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with coupling Vij = Cq/R^j [I7 2]. Including spontaneous 
emission as well as dephasing described by the stan- 
dard Lindblad operator £, the two-atom master equation 
reads 



(3) 



In a first attempt, we generalized the SARE ap- 
proach 20|, Ell to solve Eq. (j3j). But it turned out that 
the calculation of the rates in the two-atom system is 
computationally too inefficient. Also, negative rates fre- 
quently occur, prohibiting a direct Monte Carlo based 
solving technique. The latter problem can be overcome 
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FIG. 3. (Color online) (a) Rydberg excitation for a 3D sample of 500 atoms and Ce = 50000 /xm^ MHz. The hybrid model 
(dashed lines) and the single atom rate equation (solid lines) yield the same results, (b) Rydberg excitation on a ID lattice 
for different models with interaction Vnn = 2.5 MHz between adjacent atoms. The resonances at A = Vnn and A = Vnn/2 
correspond to single-atom and two-atom excitations. Other parameters are Q12 = 2 MHz, Q23 = 1 MHz, 721 = 6 MHz, 732 = 
25 kHz, r32 = r2i = 100 kHz. 



by correcting the calculated transition rates to become 
strictly positive. This eliminates the knowledge of the 
physical time evolution, but still evolves the system into 
the correct many-body stationary state [13]. 

Instead, here, we replace the time-consuming calcula- 
tion of the rates by an evolution equation based on the 



readily available steady state solution cr^'^ 
atom populations paa,bb as (see appendix) 
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By construction, it also propagates the system into the 
correct stationary state. But as it does not predict the 
physical intermediate time evolution, we drop all time 
occurrence in our equations. This allows us to iteratively 
solve for the many-body steady state with the following 
three-step Monte Carlo procedure based on the Random 
Selection Method [10]. First, a pair or a single atom is 
randomly chosen. Second, the effective detuning for the 
atom/pair is calculated based on the current configura- 
tion of ah atoms. With this detuning, the steady state 
of the atom/pair is calculated. Third, a random number 
is drawn, and the population of the atom/pair is set in a 
Monte Carlo step according to the probabilities given by 
its steady state population distribution. The three steps 
are repeated until the system converges to the many- 
body steady state. With our method, several hundred 
atoms can be simulated easily. Fitting the runtime up 
to N = 1000 results in the scaling law T ^ N^-^^ and 
a single realization including 1000 atoms takes about 0.5 
seconds on a 3.40 GHz CPU. 

In order to verify the validity of the KM we first con- 
sidered the limit of no pairs, and found results identical 
to the SARE calculations in [20]. Next, we calculated 
the Rydberg population probability pss = pill^ 
compared the data to estabhshed theoretical approaches. 
In Fig. [3l^a), we show results for 3D random samples 



of atoms with different gas densities. The suppression 
of excitation at higher densities can be well understood 
in terms of the Rydberg blockade. In all cases, our 
HM agrees perfectly with the SARE [20], indicating the 
consistency of our method. However, significant devia- 
tions occur if the geometry is changed to a ID lattice 
as shown in Fig. [IJb). Nearest neighbors interact with 
l^N = 2.5 MHz. The results in Fig.[3fb) show that next 
to the single- atom excitation peak at A = 0, the HM pre- 
dicts a resonance at A = I^n, and an additional weaker 
resonance at A = The first condition describes a 

resonant excitation process for an atom whose neighbor 
is already excited to the Rydberg state |3), because the 
laser detuning compensates the interaction energy shift 
at this point. The second and smaller resonance arises 
from a 2-photon-process |22) |33) simultaneously ex- 
citing two neighboring atoms. The SARE is not capa- 
ble of producing this second feature because it does not 
include the exact two-body interaction. To confirm the 
presence of the 2-photon-resonance we also performed the 
simulation with the cluster expansion model (CE) used 
in [14], which also shows this structure. But interest- 
ingly, all three approaches significantly deviate already at 
A = l^N- It should be noted, however, that the lattice 
geometry is a significant challenge for all models because 
of overlapping pairs and non-negligible higher-order cor- 
relations. 

We now turn to our main results on the predictive 
power of the different models for correlated many-body 
systems. For this, we calculate the pair correlation func- 
tion g^'^\r) for a disordered one-dimensional gas. It de- 
scribes the conditioned probability of having two Ryd- 
berg excitations of atoms with distance r 
We use the definition [T^- 
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(a) cluster expansion hybrid model (b) rate equation hybrid model 




detuning A [MHz] detuning A [IVIHz] detuning A [IVIHz] detuning A [IVIHz] 

FIG. 4. (Color online) Pair correlation function g^'^\r) for different laser detunings. (a) The models reproduce the Rydberg 
blockade and the uncorrelated regime. In the intermediate region resonances which indicate an excitation enhancement occur. 
The CE does not yield the line at 2 ■ Ri. The HM data was obtained from averaging 200000 Monte Carlo runs, noise is due to 
statistics, (b) The strong resonance lines are shown in detail and the process |23) |33) at Ri can be seen. The HM shows 
|22) |33) at R2 in addition. Parameters like in Fig. [3l except Ce = 900/27r /xm^ MHz,niD = 0.1 //m~^. 



Here denotes the sum over all pairs of atoms with 

distance r. For two uncorrelated atoms the pair correla- 
tion function is one. Inside the Rydberg blockade regime, 
where only one excitation can be present, it becomes zero. 

Results for g^'^\r) for different values of the laser de- 
tuning A are shown in Fig.|4j As expected, the blockaded 
and the uncorrelated regimes appear clearly for small and 
large distances in Fig. |4l^a). In between, values g^'^^ > 1 
occur for some distances and detunings. These values in- 
dicate spatial order, originating from a high probability 
for multiple excitation at selected distances. Most promi- 
nent is the resonance line at about 2 /am which can be 
characterized by the condition Ri = (Ce/A)^/^. Just as 
in the lattice simulations, it corresponds to the resonant 
excitation of a second atom when an atom with distance 
Ri is already in the Rydberg state. In this case the ef- 
fective detuning Aeff = A — Cq / Ri vanishes. Comparing 
the results for the HM and CE in Fig. Hl^a), it can be 
seen that the CE does not predict the second peak visi- 
ble around distances 4 /im and for positive detunings in 
the HM results. This peak occurs if two atoms with dis- 
tance Ri are already excited, and a third atom again with 
distance Ri to either of the two atoms is excited in addi- 
tion. The two outermost atoms in this arrangement then 
have distance 2 Ri . Note that in a higher dimensional ge- 
ometry this resonance line smears out, because then the 
distance between the first and the third atom is not nec- 
essarily 2Ri. The CE does not predict this higher-order 
spatial correlation since triply excited states are not part 
of the model. In contrast, the SARE method correctly 
predicts the correlations at 2 Ri. 

Next, we analyze the resonance originating from or- 
dering between nearest neighbors, and show a magnified 
section of Fig. HJ^a) in (b). The left panel of (b) shows 
the results obtained from the SARE, and contains the 
single resonance at Ri. But in the result of the HM, 
another resonance appears. Its condition can be deter- 



mined as R2 = [C6/(2A)]i/^ ^ 0.89 • Ri. This resonance 
again originates from the 2-photon-process |22) |33). 
Interestingly, in this case, the CE correctly predicts the 
double resonance, whereas the SARE does not. We note 
that the calculation of many-body Hamiltonian dynam- 
ics in two- level systems also shows, and thus confirms, 
the discussed resonance lines. 

In summary, we presented a hybrid model for the sim- 
ulation of large ensembles of Rydberg atoms. It is based 
on an exact two-atom calculation for pairs with distances 
below a characteristic length scale Lcr^ combined with an 
approximate treatment via effective detunings for more 
distant atoms. We proposed a method to iteratively solve 
for the steady state of the HM based on the Random 
Selection Method without the need for time consuming 
calculations of transition rates. This way, high predictive 
power is combined with fast calculation times. We found 
agreement of the hybrid approach with existing models 
for simple observables such as the Rydberg excitation 
probability. But different predictions are found for pair 
correlation function g^'^\r). The HM predicts structures 
which both, the CE and the SARE models, can only par- 
tially reproduce. We identified the additional structures 
as originating from multi-atom processes. This not only 
demonstrates the capability of the HM to characterize 
higher-order correlations, but also that higher-order cor- 
relations cannot be neglected even at low densities. 

Financial support by the Heidelberg Center for Quan- 
tum Dynamics is gratefully acknowledged. 



Appendix 

For the propagation of the system into the many-body 
steady state we use a Monte Carlo algorithm based on the 
Random Selection Method. In this supplement, technical 
details of the scheme are described in more detail and a 
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flowchart of the algorithm is shown in Fig. [5] [23|. 



Evolution equation 

Rate equation based models replace the dynamics of 
the density matrix p by effective equations for the state 
populations only. For example, for a single atom, if the 
vector containing all populations is denoted by a, then 

a = Aa , 

with a matrix A containing rates describing the transi- 
tions between the different states. For this reduction to 
the populations, the dynamics of the coherences is ne- 
glected. This in general leads to an approximate descrip- 
tion of the dynamics. But the steady state a^^^) can be 
precisely predicted, because there p^^^^ = and the dy- 
namics of the coherences is zero by deflnition. Note that 
the steady state is unique for quantum optical systems 
as considered here. 

As outlined in the main text, to avoid the problem of 
negative rates and to signiflcantly speed up our calcula- 
tion, we instead use the evolution equation (written in 
component form) 



b 

Bab = - Sab , 
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Convergence? 
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with a tailored evolution matrix B. This matrix is opti- 
mized in the sense that it can be calculated much more 
rapidly than the usual transitions rates in A, as it only 
depends on the steady state populations a^^\ At the 
same time, it preserves the total probability and evolves 
the system into the steady state, as Ba''^^^ = 0, ensur- 
ing consistency with the condition a^^^^ = Ba^^^"^ = 0. 
Eq. (4) in the main text is the straightforward generaliza- 
tion of this tailored evolution equation to the two-atom 
case. 



Monte Carlo algorithm 

1. Initialization. The positions of the N atoms are 
determined according to the desired density distribution. 
From the point of view of the algorithm, no restrictions 
apply to the positions as long as the conditions are such 
that a perturbative treatment in the order of correlation 
is meaningful. If the separation of two atoms is smaller 
than the critical distance Lcr^ they are treated as a pair. 
In total we get Np pairs and Na single atoms. Ideally one 
should have 2Np ^ Na = N. But in practice, 2 TVp + 
Na > N. For example, if three atoms are placed on 
a line with distance between neighboring atoms smaller 
than Lcr, then the three atoms count as two pairs, such 
that 2 7Vp + 7Va = 4>7V = 3. Next, the initial state of 
all N atoms is set to the ground state 

2. Choice of atom/pair. In the flrst step of the Monte 
Carlo sequence, the atom or pair is determined whose 
state will be updated in the sequence. For this, a random 
integer between 1 and Na + Np is drawn which deter- 
mines the atom or pair which is handled in the current 
Monte Carlo step. In the following, the indices labeling 
the atoms are denoted as i in case of a single atom and 
ii and 12 in case of a pair. 

3. Effective potential. Next, the effect of all other 
atoms on the atom/pair chosen in step 2 is calculated. 
For this, the interaction potential of the atom/pair with 
all other atoms is summed up and included in an effective 
detuning. For a single atom i it reads A^^ = A — Vij 
with the sum running over all other atoms j which 
are in the Rydberg state |3). In case of a pair two ef- 



fective detunings A 



eff 



ii2) 



Y^'j-Li^ Vi^j have to be calculated. 



E'j^i, Viij and A^^ 



Note that the 

interaction Vi^i^ is included exactly in the two- atom de- 
scription of the pair chosen in step 2. 

4. Steady state. Next, the steady state of the 
atom/pair chosen in step 2 under the action of the ef- 
fective potential calculated in step 3 is determined. For 
the single atom the steady state for the population 



Pa 



(Saa) is denoted by aa^^\ For a pair of atoms ii 



JSS) 



and Z2, the steady state of Paa,bb = (^da is labeled 

as cr^5^^ Here, a,6G{l,2,3}. This step is the most time 
consuming part of the simulation and should therefore 



FIG. 5. Schematic algorithm for the hybrid model. 
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be highly optimized. By comparing different approaches, 
we found the following method to be the fastest numeri- 
cally stable approach. The master equation is linear and 
can be written as p = M • p where M is a sparse ma- 
trix of the size 9x9 or 81x81 for single atoms or pairs, 
respectively. The fact that the master equation is Her- 
mitian can be exploited to obtain a real matrix M. The 
problem of finding the steady state is equivalent to find 
the eigenvector of M with eigenvalue 0. Since all states 
are coupled, the dimension of the corresponding eigen- 
system is one and the steady state is unique. M is now 
decomposed into a product of an orthogonal and an up- 
per triangular matrix in a QR decomposition [2J]. Since 
the matrix M is sparse, Givens Rotations are the best 
approach to transform M successively into an upper tri- 
angular matrix R [24] . R does not have full rank and the 
last element on the diagonal vanishes. The steady state 
can then be obtained by back substitution. In practise it 
is useful to arrange the elements of the vector p such that 
the populations are located at the end, because then the 
back substitution can be stopped earlier. 

5. State update. Next, the state of the atom/pair de- 
termined in step 2 is modified according to the population 
probabilities determined in step 4. For this, a random 
real number r between and 1 is drawn. In the case 
of a single atom, the state is changed to |/), where / is 
the largest integer with Xl/c=i < ^- ^^^r pairs, the 
procedure is analogous. Note that the new state can be 
identical to the old state. 

6. Loop of the Monte Carlo sequence. The steps 2-5 
are now repeated until the system is converged. We have 
found that for typical parameters, approximately 10 • 
steps are sufficient to ensure convergence. 

7. Calculation of observables. Finally, the observables 
can be evaluated. Typically, an average over many spa- 
tial realizations and Monte Carlo trajectories is neces- 
sary to obtain good statistics and smooth results ex- 
pected from larger ensembles of Rydberg atoms. Note 



that throughout the Monte Carlo evolution, each atom 
is in a definite atomic state at all times, and the diagonal 
elements of the density matrix only originate from the 
ensemble averaging. In our work we focus on the evalu- 
ation of g^'^\r)^ and Eq. (5) in the main text shows how 
to calculate it from the simulation results. 
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